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We investigate the influence of curvature and topology on crystalline wrinkling patterns in generic 
elastic bilayers. Our numerical analysis predicts that the total number of defects created by adiabatic 
compression exhibits universal quadratic scaling for spherical, ellipsoidal and toroidal surfaces over a 
wide range of system sizes. However, both the localization of individual defects and the orientation 
of defect chains depend strongly on the local Gaussian curvature and its gradients across a surface. 

Our results imply that curvature and topology can be utilized to pattern defects in elastic materials, 
thus promising improved control over hierarchical bending, buckling or folding processes. Generally, 
this study suggests that bilayer systems provide an inexpensive yet valuable experimental test-bed 
for exploring the effects of geometrically induced forces on assemblies of topological charges. 


Topological defects and geometric frustration are in¬ 
trinsic to two-dimensional (2D) curved crystals [T]. The 
minimal number of defects in a periodic polygonal tiling 
is dictated by Euler’s theorem [5], which relates the sur¬ 
face geometry to the total topological defect charge. The 
hexagonal soccer-ball tiling is a canonical example, re¬ 
quiring 12 pentagonal defects that are also realized in 
Ceo fullerenes [5]. However, large 2D crystals often ex¬ 
hibit defect numbers that go substantially beyond the 
minimal topological requirement |4]. These excess de¬ 
fects aggregate in molecule-like chains m that relieve 
elastic energy costs arising from a mismatch between the 
crystal symmetry and the curvature of the underlying 
manifold B ini- The aggregation of curvature-induced 
defects plays an essential role in various physical pro¬ 
cesses, from the classic Thomson problem of distributing 
discrete electric charges onto a sphere m to the assem¬ 
bly of virus capsules m and the fabrication of colloido- 
somes [T^ , toroidal carbon nanotubes m, and spherical 
fullerenes from graphene mi. Over the past decade, con¬ 
siderable progress has been made in understanding crys¬ 
tal formation in spherical 00115 ] and more complex ge¬ 
ometries [TBfl22| . Yet, empirical tests of theoretical con¬ 
cepts have remained restricted Emin to paraboloids 
or mean-curvature surfaces, owing to the lack of tractable 
experimental model systems. 

Here, we show through theory and simulations that 
curved elastic bilayer materials offer a promising test-bed 
for studying defect crystallography in arbitrarily shaped 
2D geometries. Building on a recently derived and ex¬ 
perimentally validated scalar field theory mi, we first 
confirm that thin-film wrinkling reproduces previously 
established results for the crystal formation on spheri¬ 
cal surfaces 0[26]. Subsequently, we demonstrate how 
curvature and topology determine elastic defect localiza¬ 
tion on surfaces with non-constant curvature. For typi¬ 


cal experimental parameters |26| , our analysis reveals the 
emergence of previously unrecognized robust superstruc¬ 
tures, suggesting the usage of topology and geometry to 
control defect aggregation. 

Our elastic bilayer system comprises a thin stiff film 
adhered to a soft curved substrate. Recent experi¬ 
ments [26U28] with spherical substrates showed that, un¬ 
der weak compression, such films can wrinkle into a crys¬ 
talline dimpled pattern (Fig. [^. The experimental pat¬ 
terns are described quantitatively by a generalized Swift- 
Hohenberg theory [5^, which is employed here to ob¬ 
tain predictions for more general geometries. Measur¬ 
ing lengths in units of film thickness h and focussing on 
leading-order effects, the theory relates the normal dis¬ 
placement field u of the deformed film to the minimum 



FIG. 1. (color online) Stress-induced crystalline wrinkling 
patterns in a thin film of thickness h adhering to a soft core 
(top), obtained by minimization of Eq. Q, and their dual 
hexagonal Voronoi tessellations (bottom) for different surface 
geometries: (a) sphere {R/h — 70), (b) ellipsoid {Rx = 2Ry = 
2Rz = IlOh), (c) torus {r/h — 16, R/h = 80). The color bar 
represents the surface elevation. The outlined surface domain 
in (b) highlights chains of defects. Voronoi cells are color- 
coded by their coordination number Z. 
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FIG. 2. (color online) (a) The total number of defects grows linearly with y/N, where N is total nnmber of lattice units, 
exhibiting similar slopes for all geometries (solid line: linear ht for spheres). Inset: Excess dislocations on spheres were 
predicted [6] to increase linearly with reduced radius p = R/X for colloidal crystals (solid line). For comparison, the best-fit 
power law (p — pc)^ with pc = 4.5 ± 0.4 and /3 = 0.67 ± 0.08 is also shown, (b-d) The total defect charge grows differently with 
integrated Gaussian curvature I for different geometries. The gray-shading represents the conditional PDFs of the total charge 
Q for a given value of I. The red dashed line corresponds to linear growth I = [tt/Z)Q. Dark regions in the surface sketches 
illustrate integration domains. 
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where k = Ef/{1 — i/^), and Ef is the Young’s 
modulus of the film with undeformed surface element 
doj. The Poisson’s ratio ly is assumed to be equal for 
the film and substrate. The nonlinear term r(u) = 
[(I — v)b°‘^VauVpu + 2v'H{VuY\ u represents stretch¬ 
ing forces, with surface gradient V and Laplace-Beltrami 
operator A. The trace "H of the curvature tensor 6 “^ 
dehnes the mean curvature and JC denotes the Gaus¬ 
sian curvature. Surface wrinkles form when the film 
stress exceeds the critical buckling stress, corresponding 
to a < Oc = 37 q where 70 dehnes the pattern wave¬ 
length A = 27r/-\/6|7o|, while the amplitude is controlled 
by c [ 55 . The curvature-dependent r(u)-term deter¬ 
mines the symmetry of the predicted wrinkling patterns. 
In the planar case with T{u) = 0, Eq. ([^ reduces to the 
Swift-Hohenberg model [521 minimization of Eq. ([^ 
produces stripe-patterns. Eor r(M) yf 0, the symmetry 
u —> —u is broken, causing a transition to hexagonal 
dimple-patterns |5S]. A systematic derivation [5S] spec- 
ihed all parameters in Eq. ([^ in terms of known mate¬ 
rial and geometric quantities, so that our predictions can 
be directly compared to future experiments. Our simu¬ 
lations use the material parameters of Ref. [25] for the 
hexagonal phase with A ^ 9.1 throughout. 

We analyze Eq. Q for: spheres of radius i?, prolate 
ellipsoids of principal axis Rx = 2Ry = 2Rz, and tori 
with aspect ratio r/R between major axis R and minor 
axis r (Fig.[^. Eor each geometry, we choose a/oc « 0.98 
and add small random perturbations to the unwrinkled 
him uq = 0. We then minimize Eq. 0 numerically us¬ 
ing a custom-made hnite element method |30] . Erom the 


stationary displacement held, the dimple centroids are 
identihed as lattice points and the corresponding Voronoi 
tessellations are constructed (Fig. [^. We dehne the 
topological charge s = 6 — ^ for each lattice point, where 
Z is the coordination number. Probability density func¬ 
tions (PDFs) and statistical averages are obtained from 
simulations with different random initial conditions but 
otherwise identical parameters. 

Since the total number of lattice units N is propor¬ 
tional to the surface area, '/N is a geometry-independent 
measure of the system size. We now characterize the to¬ 
tal number of defects and hnd that it grows linearly 
with for all geometries (Fig. [^), while generally 
exceeding the Euler bound for the minimal number of 
defects. For example, we observe exactly 12 topologi¬ 
cally required defects of charge s = -1-1 {positive discli- 
nations) only for small spheres, whereas increasing 
linearly with slope m = 4.5 ± 0.6 above the critical size 
yAc = 14.2 ± 3.6 (Fig. 0 ) . In terms of the reduced ra¬ 
dius p — R/X, Nc corresponds to a critical value p* ^ 4.3, 
in good agreement with experiments |26j . 

The increase of with system size can be explained 
as follows. Each disclination imposes a set change of 
Gaussian curvature, independently of -s/iV. If the mis¬ 
match with the substrate’s target curvature becomes too 
large, additional defects are introduced to screen curva¬ 
ture, thereby lowering the total elastic energy PEI]. To 
preserve the total charge, excess defects appear as neu¬ 
tral pairs of opposite charge called dislocations. For large 
systems, defects typically form longer chains classified as 
either neutral pleats or charged scars (Sj. For spheres, 
the number of defects per scar is predicted to grow lin¬ 
early above pc with slope « 0.41p PED- This scaling has 
been experimentally verified in colloidal crystals P , and 
agrees with our wrinkling simulations, although a power- 
law ~ (p — Pc)^ with Pc = 4.5 ± 0.4 and {3 = 0.67 ± 0.08 
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also fits our data well (Fig. inset) |32]- These find¬ 
ings illustrate that geometry-induced defect formation is 
insensitive to the details of the lattice interactions [33], 
corroborating that wrinkling patterns provide a viable 
model to study generic aspects of curved crystals. 

The screening effect of dislocations and its dependence 
on geometry become manifest in the relationship be¬ 
tween Gaussian curvature and topological charge 0 HH] ■ 
The Euler and Gauss-Bonnet theorems connect the sum 
of topological charges Si for all elements, Q = 
with the surface integral of the Gaussian curvature, 
I = JAlCdA = 7r/3 Si = 27rx, where x is the Euler 
characteristic of the surface. How well this relationship 
is satisfied over a subregion of the surface provides in¬ 
sight into the geometry dependence of defect localization. 
For spheres, our results are consistent with Gauss-Bonnet 
theorem; Q increases linearly with I (Fig. [2^)- By con¬ 
trast, for ellipsoidal geometries, Q grows faster than / 
near the poles (|0| = tt/2) and there is an accumulation 
of positive charges in high-curvature regions (Fig. i=)- 
Although tori require no topological charge (x = 0), our 
simulations predict the creation of defects that help the 
dimpled crystal comply with the curved substrate geome¬ 
try M- In the outer region of the torus, where Gaussian 
curvature is positive, we find that Q grows faster than 
linearly with 7, which is qualitatively similar to the el¬ 
lipsoidal case but with larger spread. 
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FIG. 3. (color online) Curvature-induced defect localization 
on ellipsoids, (a) Isolated pentagonal -|-1-defects accumulate 
in high curvature regions; p is the average number of single 
defects per ellipsoid, (b) Although defect chains form prefer¬ 
entially near the equator (|(^| = 0) their orientation angles a, 
measured relative to the tangent vectors t, show no signifi¬ 
cant ordering. (c,d) Voronoi tessellations for Rx/h = 40 and 
Rx/h — 160. Ellipsoids with Rx/h > 110 show a weak align¬ 
ment of the lattices along lines = 0. 


critical transition angle 


Another striking phenomenon is the curvature-induced 
localization and segregation of oppositely charged de¬ 
fects. For ellipsoids, we find that the PDF for angu¬ 
lar position of positively charged disclinations increases 
strongly towards the poles, where the local Gaussian cur¬ 
vature is large enough to support them (Fig. [^). With 
increasing size, additional scars and pleats appear. Their 
centroid positions cluster in the equator region around 
(j) = 0 (Fig.j^), where the curvature is low and thus can¬ 
not support isolated disclinations. To study the orien¬ 
tations of these extended defect ‘molecules’, we measure 
the orientation angle a enclosed by the end-to-end vec¬ 
tor V and the tangent t along lines with Q = const. We 
find no significant orientational order for positive or neu¬ 
tral chains (Fig. -d), consistent with earlier simulations 
based on an inflation packing algorithm [34j . 

Tori, in contrast to ellipsoids, contain regions of pos¬ 
itive and negative curvature and are more prone to 
striped wrinkling. To identify the conditions for the pure 
crystal phase, we recall that hexagonal patterns require 
r(M) ^ 0 in Eq. 0 whereas local stripe solutions emerge 
for r(u) —>•0 [ 23 . The phase boundary can be estimated 
by first parametrizing the torus using the standard coor¬ 
dinates {(f>, 9), then assuming a stripe-like wrinkle pattern 
symmetric along 6 and finally inserting u{9,(f)) = 
into the condition r(w) = 0. Solving for (j)^ we find the 
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which is independent of the system size and holds as long 
as 7?, r ^ A. For rubber-like materials (iz « 0.5) solutions 
±(j)c G [f,7r] exist for aspect ratios r/R > 1/3. We thus 
expect stripe-like wrinkles to dominate near the inner 
rims of thick tori. For r/R < 1/3, (pc is imaginary; the 
symmetry-breaking r(u)-term is globally non-zero, im¬ 
plying a purely hexagonal phase throughout the torus. 
To verify these predictions, we performed simulations for 
tori with 0.2 < r/R < 0.8. Defining (j>c in simulations 
as the angle beyond which less than half of the wrinkled 
surface takes the form of hexagonal dimples, we find good 
agreement with Eq. see Eig. |^. The existence of a 
pure hexagonal phase for r/R <1/3 provides a design 
guideline to study toroidal crystals in future wrinkling 
experiments. 

Eocusing on crystalline patterns on slender tori with 
r/R = 0.2, our simulations show that defect localiza¬ 
tion is strongly controlled by the interplay of Gaussian 
curvature and topology. The requirement of a vanishing 
net charge implies that positive and negative disclina¬ 
tions appear in pairs. Their spatial arrangement can be 
rationalized with an electrostatic analogy |18j . in which 
defects are interpreted as charged particles and curva¬ 
ture acts as an electric field. In this picture, five-fold 
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FIG. 4. (color online) Defect localization and superstructures on tori, (a) Pure hexagonal wrinkling phases are stable only for 
tori with r jR < 1/3 (shaded). Stripe-like patterns occur at angles 4> > 4>c for larger aspect ratios. Numerical simulations agree 
well with the theoretical predictions (Eq. 2, solid line), (b) Isolated penta- and heptagonal defects segregate; p and n are the 
average numbers of positive and negative disclinations per torus, (c) Orientations of charged scars and neutral pleats correlate 
strongly with their transversal centroid position \<j)\. (d) Defect positions are also strongly correlated along 6. (e,f) While this 
positional ordering is less prominent in small systems [(e), R/h = 40], it becomes apparent for larger systems [(f) R/h — 160]. 
The toroidal superstructure (solid black lines) appears to follow geodesics of minimal integrated Gaussian curvature (dashed 
green line). 


disclinations are attracted to regions of positive Gaus¬ 
sian curvature at the outer rim of the torus (</> = 0), 
whereas seven-fold disclinations migrate to the inner re¬ 
gion (|(/)| = tt). This geometry-induced separation of 
charges is directly reflected in the PDF of the individ¬ 
ual disclinations (Fig. &)• Analogously to the ellip¬ 
soidal case, the total number of isolated disclinations 
(see the average numbers p and n for positive and neg¬ 
ative charges in Fig. |^) decreases with system size, as 
defects tend to aggregate in chains (Fig. |^,f). Inter¬ 
estingly, we find that the electrostatic analogy extends 
to defect chains: positive scars screen Gaussian curva¬ 
ture on the outer rim; negative scars appear in the in¬ 
ner region; and neutral pleats concentrate in the regions 
of vanishing Gaussian curvature, \(j)\ = 7r/2 (Fig. &)■ 
Defect chains also become oriented by geometric forces. 
Measuring the orientation angle a of a chain relative to 
the tangent vector t along the ^—direction (Fig. 
we find that charged scars preferentially align parallel 
to the equatorial lines such that a ^ 7r/2, whereas neu¬ 
tral pleats tend to orient vertically with a ~ 0 (Fig. &)■ 
This ordering can be understood qualitatively by con¬ 
sidering the end-points of a defect chain. For scars, both 
end-points have the same charge and therefore migrate to 
regions of same Gaussian curvature ((/> ^ 0 for positively 
charged endpoints, ~ tt for negatively charged ones), 
effectively orienting the scar perpendicular (a = 7r/2) 
to lines (j) = const. By contrast, pleats have oppositely 
charged ends and hence mimic electric dipoles that be¬ 
come oriented by curvature to achieve a = 0. 


Remarkably, our simulations reveal that defects not 
only orient and segregate in the torus curvature field - 
they also break the rotational symmetry of the toroidal 
crystal in favor of an emergent discrete symmetry. More 
precisely, by analyzing the lattice structure on large tori, 
we find that defects arrange along an undulating periodic 
deformation pattern (highlighted through black lines in 
Fig. 1^). This superstructure still carries a fingerprint 
of the underlying toroidal geometry: geodesics = 
{(j), 9{(j))) on a torus are solutions of [35] 


dfj} {R + r cos (j))^y{R + rcoscj))^ — 

where c is a constant obeying Claireaut’s geodesic rela¬ 
tion c = (i? -k r cos (/)) sin'0, with ip denoting the angle 
between the tangent g' and t [35]. Integrating Eq. ^ 
numerically, we find that the lattice deformation follows 
a geodesic passing through p = 1.68 at its highest point 
with ip = TT 12 (dashed green line in Fig. 1^). The phase 
6q, found by matching the phase of the geodesic with 
that of the lattice, varies between samples. The inte¬ 
grated absolute Gaussian curvature along this curve is 
minimal among all geodesics that do not wind around the 
(/)-coordinate. We therefore hypothesize that alignment 
along a geodesic of minimal absolute Gaussian curvature 
yields an energetically favorable crystal structure, for at 
least those lattice parts close to where the geodesic re¬ 
main nearly straight - much like wrapping a torus with 
an inextensible ribbon. If this geometric argument is 
correct, the superstructure should become independent 
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of the lattice size, as, R/h increases. To test this hypoth¬ 
esis, we simulated system sizes R/h G {80,120,160}, and 
measured the angular distance A0 between positive scars 
near the outer rim. These simulations indeed reveal a 
peak in the PDF P{A6) near 7r/4, which corresponds to 
one quarter of the geodesic wavelength, independently 
of system size (Fig. |^, black curves). Moreover, nega¬ 
tive scars on the inside of the torus appear in phase with 
positive scars {A6 = 0), while neutral pleats arrange be¬ 
tween scars (A0 ~ tJ'/S). Thus, the geometric lattice 
superstructure also controls defect-chain localization. 

A significant advantage of elastic crystals over fluid- 
based systems is their fabrication versatility, which en¬ 
ables the exploration of arbitrary geometries and topolo¬ 
gies. Our results show that spatially varying curvature 
can lead to emergent superstructures that determine de¬ 
fect localization. Since defects trigger secondary instabil¬ 
ities |36j . this previously unrecognized phenomenon can 
be exploited to control hierarchical buckling and folding. 
More generally, our analysis implies that elastic crystals 
provide a rich model system for studying the profound in¬ 
terplay between geometric forces and topological charges. 
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